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Abstract 

We revisit the coherent or delocaHzed to self-trapping transition in an interacting bosonic quan- 
tum fluid confined in a double well potential, in the context of full quantum calculations. We 
show that an AT-particle Bose-Hubbard fluid reaches an stationary state through the two-body 
interactions. These stationary states are either delocalized or self-trapped in one of the wells, the 
former appearing as coherent oscillations in the mean-field approximation. By studying one- and 
two-body properties in the energy eigenstates and in a set of coherent states, we show that the 
delocalized to self-trapped transition occurs as a function of the energy of the fluid, provided the 
interparticle interaction is above a critical or threshold value. We argue that this is a type of 
symmetry-breaking continuous phase transition. 



I. INTRODUCTION. 



The experimental realization of a single bosonic Josephson junction in two weakly linked 
Bose- Einstein condensates (BEC) has revealed fundamental features of the dynamical nature 
of an ultracold interacting Bose gas in a double-well potential[T]. The atomic transport has 
shown the existence of Josephson tunneling or macroscopic quantum self-trapping regimes 
depending on the interatomic interactions and on the initial populations in the wells. These 
opposite states are similar to the well known superfiuid and Mott insulator states observed 
in optical lattices [21 E] in the sense that the quantum transition among them is the result 
of the non-linear condensate self-trapping interactions. 

The theoretical analysis of this system has been based on the Bose-Hubbard Hamiltonian 
in the two-mode approximation [H |5l El [3 [H [9] . Besides providing a fairly good description of 
the experimental situation [Ij, this model has been extensively studied on its own, specially 
within a mean-field or semiclassical approximation. This model has lead to significant un- 
derstanding of the richness of the physical problem at hand; additionally, it can incorporate 
asymmetric two-well potentials and external time-dependent driving fieldsjl]. 

Beyond the studies within the mean-field approximation, this model is amenable to nu- 
merical full quantum calculations, so far up to = 1000 particlesjH IH [U |9]. These results 
and mean-field calculations agree for a short time while for longer times the full quantum 
approach shows that the system reaches a state, that one may call "statistically stationary" 
as we shall specify below, punctuated by coherent revivals. Here, we revisit the symmet- 
ric two-mode Bose-Hubbard model within numerically exact quantum calculations to, first, 
verify the previous results and, second, to accomplish a more detailed description of the self 
trapping transition. 

The mean-field approach, which approximates the A^-body dynamics by a set of non-linear 
coupled dynamical equations, is limited and incomplete since it describes the evolution of the 
expectation values of one-body operators only. However, from the full quantum dynamics 
one may, in principle, enquire about the dynamics of one-, two-, three-, and up to A^-body 
properties. Due to the relative simplicity of the present model one can calculate either the 
full A^-body wave function for a pure state or the full A^-body density matrix for a mixed 
state (for a system up to A^ = 10,000 particles, say). Thus, full information of the state of 
the system may be obtained. 
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Anyhow, from a macroscopic point of view, most of the measurable thermodynamic 
and/or many-body properties, such as energy, temperature and Green's functions in general, 
are given in terms of one- and two-body operators [TU|. These in turn are exhaustively 
described by the one- and two-body reduced density matrices. Although we are able to 
calculate the complete one- and two-body density matrices from the exact solution, here for 
purposes of showing our main points we shall limit our study to the following measurable 
quantities: the number of particles in one of the wells, given by the expectation value of 
a one-body operator Ni, and a tunneling correlation function, the expectation value of a 
two-body operator J^Jx- Below we give the explicit form of these operators. 

Our results may be summarized as follows. Since we are dealing with a macroscopic 
system, following the principles of statistical mechanics, we choose to describe the state of 
the system by specifying the number of particles and the expectation value of the energy e; 
that is, we shall work in a microcanonical ensemble. For this to be meaningful, we must limit 
ourselves to states whose energy mean-square deviation is small; in other words, to states 
localized in energy but not necessarily the energy eigenstates. For purposes of calculations 
we use a well-known family of coherent states [HI [12] such that the previous requirement is 
satisfied. We point out here that the particular state with all particles initially in one of the 
wells, analyzed in most of the studies with full quantum solutions|ll El |9], belongs to the 
family of such coherent states. 

The first important conclusion is that, for all the above mentioned initial states and for 
all values of the strength of the particle interactions, the system reaches a state that can be 
called "statistically stationary". That is, the distributions of the time- dependent expectaiion 
values < iVi > and < JxJx > are peaked around a value dependent only on A^ and e, with a 
width that becomes relatively smaller as A^ is increased; in other words, the distribution of 
those values is sharply peaked. Secondly, we find that the mean values of those distributions 
follow the same function as the expectation values of the operators A^i and JxJx in the 
energy eigenstates, truly stationary states. 

For small values of the interaction strength, typically, the values of < iVi > and < JxJx > 
oscillate following the coherent oscillations predicted by the mean-field calculation for a 
brief period of time, then "decay" to a stationary value for a long period of time (e.g. 
< A^i >^ N/2) followed by revivals. As the number of particles is increased, the time spent 
in the revival regions and the time spent in the stationary region both increase. However, 
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the ratio of the latter to the former also increases, such that most of the time is spent in the 
stationary region. Could this process be called decoherence? It certainly so if experimentally 
one is able to measure few (i.e. one- and two-) body properties only. After all, this is what 
statistical physics asserts: any isolated system left undisturbed for a long time relaxes to 
a stationary (equilibrium) state. For larger values of the interaction not only the decay to 
stationary states occurs, but also the transition to self-trapped states appears. The main 
result is that there exists a critical value of the interaction, such that if the interaction is 
greater than this value, the self-trapping transition occurs as a function of the energy of the 
system. We recall that since the one-particle spectrum is finite (i.e. just two states), the 
iV-body set of allowed energy values is bounded from above and below. 

This article is organized as follows. We briefly present the Hamiltonian of the system 
and the quantities to be analyzed. We then show results that display the existence of the 
statistically stationary states and we discuss the transition from delocalized to the self- 
trapped states. We conclude with a discussion of our results and additional aspects of the 
phenomenology of this system, such as its characterization in different ensembles. 

II. THE MODEL, OBSERVABLES AND EIGENSTATES PROPERTIES 

We model the interacting quantum fluid confined in the double well through the Bose- 
Hubbard Hamiltonian considering the two- mode approximation [5] 




where A is the tunneling frequency of the two lowest energy modes in a symmetrical two 
well potential and U = Anh'^a/m represents the effective particle-particle interaction strength 
written in terms of the (positive) s-wave scattering length a. We shall specify the different 
regimes by the adimensional interaction parameter A = NU/hA. The model may also be 
seen as the collection of atoms with two internal states in an external static field, and 
whose pair interaction is long range and only acting if the particles are in the same internal 
state. We shall return to this version in the last section. We use units with h = A = m = 1. 

As it has been already established, for particular initial states the system described by 
Hamiltonian ([T| exhibits a continuous transition from particle coherent oscillations to a self- 
trapping regime when the parameter A is increased from zero to above a critical value Ac. 
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This transition may also occur for a fixed value of A, by varying the initial state, which 
is equivalent to varying the energy of the system. As a matter of fact, the experiment by 
Albiez et al.pQ using Bose atoms at very low temperatures, reports the observation of a 
self-trapped state as a function of the initial population imbalance. 
It proves convenient to introduce the one-body operators|5], 

Jx = \ (blh + 4&i) 

I = l{blbi-blb2) (2) 

that obey the SU(2) commutation relations, [Jk, Ji] = iSkimJm, and the number of particles 
N = b[bi + b\b2- With these identifications the Heisenberg equations of motion for the Jk 
operators may be transformed to their mean field counterparts by identifying —>■ NK^ 
with Kk a c-number function of time. It is believed, that such an approximation is valid for 

— s> oo, f/ — s> 0, but A finite |4j. As pointed out in the introduction, this procedure gives 
information, within the limit considered, on one-body properties only. Nevertheless, this 
approximation is susceptible of extensions to systems with external time-dependent fields, 
whose implementation in exact (numerical) quantum calculations may be quite difficult. 

By solving for the full A^-body wavefunction \ipNit)) one can calculate the A^-body density 
matrix p^^\t) = \ijjiy(t)){'ipj\f{t)\, and with this, the reduced one- and two-body density 
matrices, p^^^ and p^'^\ whose matrix elements are, 

pS = Tr {K p("^ bl) , (3) 

and 

pi? = Tr {b,K p(^) , (4) 

with the Greek subindices taking the values 1 and 2, and a = (pu) and b = {nf]), that is, the 
values (11), (12) and (22). Knowledge of these two operators suffices to determine all one- 
and two-body properties of the system. Although, we can calculate all the matrix elements 
of p^^^ and p^^^ as a function of time, we find more illustrative to limit ourselves to two 
properties, one is the number of particles in well 1, namely, 

iVi = b\bi 

= l^ + J^^ (5) 
5 




FIG. 1: (Color online) Eigenstate expectation value of the particle population in well 1 {Ni)/N as 
a function of the energy eigenvalues for A = 0.1, A = 1.0 and A = 10. Green, blue and red dots, 
correspond to = 100, 1000 and 10000 particles, respectively. 

a one-body operator; and the other is a "tunnehng correlation", 

C = JxJx, (6) 

a two-body operator. 

Before studying the time evolution of initial coherent states, we analyze the properties of 
the eigenstates of the system. First, we numerically solve the eigensystem, 7t^(A)|0„(A)) = 
£:„(A)|0„(A)), where we have indicated that we have a different set for a given A. For 
particles there are A^ -|- 1 eigenstates, the energy of the system is finite and bounded by the 
lowest eo{A) and largest eiv(A) eigenenergies. We exhibit all our results, always, for given 
values of those parameters and as a function of the expectation value of the energy in the 
given state. In Figs. 1 and 2 we show the values of < iV'i(A) >„= (0r4(A)|iVi|0„(A)) and of 
< C(A) >n= (0n(A)|C|0„(A)) as a function of the energy (0 (A)|7^(A)|0„(A)), which in this 
case equals e„(A). 

/^From Figs. [T] and [2] one can appreciate the intrinsic nature of the stationary eigenstates 
as A is varied. For A = 0.1 the system does not have any transition, while for A = 1.0 
and A = 10.0, the appeareance of the transition to the self-trapped state is evident. We 
note that this is a transition as a function of the energy for given A, that is, there is 
an energy Ec above which the system self-traps. In Fig. [T] the transition is seen as a 
"symmetry breaking" of the eigenstates since the population becomes unbalanced, while 
in Fig. |2] the two-body correlation shows a peak distinguishing between "delocalized" and 
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FIG. 2: (Color online) Eigenstate expectation value of the two particle correlation C/N'^ = 
{JxJx) /N"^ , as a function of the energy eigenvalues En for A = 0.1, A = 1.0 and A = 10. Green, 
blue and red dots, correspond to N = 100, 1000 and 10000 particles, respectively. 

self-trapped states. For = 1000 and N = 10000 particles the result is essentially the 
same. As described in previous studies [HI E], in general, when tunneling and interaction 
energies become comparable NU ~ hA, the system is prone to suffer the transition. Here, 
we find that the transition is continuous, and that occurs at a "critical" value A^, ^ 0.6 for 
A^ = 1000 particles. All this is best summarized in the phase diagram shown in Fig. |3} As 
we shall show in the next section, for A > Ac there exists a "critical" energy edA) such that, 
if the energy of the system is above it, e > Sc, the system is in the self-trap region, while if 
its energy is below the transition the system will decay to a delocalized state, regardless of 
its initial population imbalance. We have dubbed "delocalized" states as opposite to self- 
trapped. In the mean-field approach these states are called "coherent oscillations" , however, 
as described in the next section, the full quantum solution shows a "decay" to stationary 
states with half of the population in each well, thus justifying its naming. 

III. STATISTICALLY STATIONARY STATES 

As already pointed out in the literature [H El E], the dynamical study of the interacting 
Bose system with Hamiltonian ([T]) reveals, in addition to the opposite coherent and self- 
trapped states, stationary states with recurrent revivals which mean-field treatments are 
unable to capture. The existence of these stationary states is due, both, to the pair-particle 
interaction energy and to the fact that we are looking at expectation values of few-body 
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FIG. 3: Phase diagram e/N vs A for N = 1000 particles. The allowed states are those within the 
blue lines, eo(A) and £n{K). The red line is the transition from delocalized to self-trapped states, 
ec(A) for A > Ac « 0.6. The states with energy greater than ec(A) are self-trapped. 

operators. The revivals are consequence of the finite number of eigenstates comprising any 
conceivable state in this system. To exhibit this phenomenolgy we consider the following 
family of coherent states pT| [T2]. 

(\ 1/2 
sin^-"i {e/2) cos"i {9/2) e-^^^-^^^-^Ini, - ni) (7) 
Til j 

where rii and N — rii are the number of particles in wells 1 and 2, and |ni, — ni) are the 
(A^ -|- 1) Fock states. The coherent states are typically sharply localized both in the Fock 
states, as well as in the eigenenergy basis. 

For given values of the number of particles A^, we characterize the coherent states by 
their energy expectation value, e{6, (p) = {9, 0| H \9, cj)). Clearly, one can have the same value 
of e{9, 0) for different values of 9 and 0. Thus, to simplify the analysis, we choose = and 
vary 9. As mentioned in the Introduction, the state with the A^ particles in well 1 is given 
by = TT and = 0. Taken any of these states as the initial one for the system, we evolve it 
numerically, \9,(j)]t) = exp{—iT-lt / h)\9 , (p) and calculate the expectation values of Ni and C, 

N,{t) = {9,<j)-t\N,\9,<P;t) (8) 
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and 

c{t) = {e,<p-t\jj^\e,<p-t). (9) 

Figs. |4] and [5] show the time evolution of the expectation values Ni{t) and C{t) in the 
coherent state 9 = tt and = 0, for A'^i(O) = A^. Figs. [3^ and|4^ are for A = 0.1. The system 
shows coherent oscillations initially, then it falls into a stationary value (with Ni{t)/N = 0.5) 
to have revivals at almost periodic intervals. We have numerically verified that the time spent 
in the stationary regions becomes longer than the time spent in the coherent-oscillation 
sections, as the number of particles increases. Figs. |3)d and |4)d correspond to A = 1.0, 
which are within the crossover transition to the self-trapped states. There are no longer 
clear coherent oscillations and the mean value of Ni/N is already above 0.5. In Figs. ^ and 

corresponding to A = 10.0, the system is well into the self-trap region and the values of 
Ni/N and C/N"^ are essentially constant. Below, we explore the whole energy region for the 
same values of A. The purpose of these figures is to exemplify two important conclusions of 
our analysis. First, based on the facts that macroscopically and experimentally we are very 
much limited to access physical properties of few bodies only, such as energies, temperatures 
and correlations functions, then, if we can make measurements of, say, A^^i and C at different 
times (starting each realization in the same state), one necessarily arrives to the conclusion 
that the system is, within fluctuations, around a mean stationary value. We shall call 
these state statistically stationary states to distinguish them from the true stationary energy 
eigenstates. This is clearly in agreement with the general principles of statistical physics. 
We have verified that the behavior shown here for the coherent state 6 = it and = 0, is 
typical for any value of 6 and 0. The second aspect we want to allude is the transition from 
coherent oscillations or delocalized to self-trapped states as the pair-particle interaction is 
increased. This point, however, is much better understood with the analysis we provide 
below. 

Once we have shown that the system reaches statistically stationary states, we can de- 
scribe its properties by the given values of the mean of A'^i(t) and C(t), as well as their 
standard deviation, for a whole family of coherent states. This family is built in the fol- 
lowing manner. For given values of and A, and taking = 0, we find the values of 6 
such that the expectation values of the energy in the coherent states span the whole interval 
£^o(A) < e(6',0) < eAr(A). Clearly, the value of 9 for a given energy is not unique. We shall 
see that, depending on this value, the coherent state can break the symmetry. That is, it 

9 



A 0.6- 



A = 0.1 



20000 



_i_ 



_i_ 



40000 60000 
t 



A = 1.0 



A = 70.0 




20000 40000 



60000 



20000 



40000 



60000 



FIG. 4: Time evolution of the expectation value of the number of particles in well 1, Ni{t)/N, in 
the coherent state 6* = tt and </) = 0, for iV = 1000 and A = 0.1, 1.0 and 10.0. 
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FIG. 5: Time evolution of the expectation value of the tunneling correlation, C{t)/N'^, in the 
coherent state 9 = ir and (p = 0, for N = 1000 and A = 0.1, 1.0 and 10.0. Note the scale in y-axis 
for cases A = 0.1 and 10.0. 

localizes either near the well 1 or the well 2, once the condition for self-trapping is satisfied, 
i.e. for A > A^. 

Figs. [6] and [T] show the mean values of Ni{t) and C{t) for = 1000 and for A = 0.1, 1.0 
and 10.0 as a function of the expectation value of the energy of the state. In the same graph 
we have included the eigenstates expectation values of Figs. [T]and[2| In general, we see that 
the mean values from the coherent states agree fairly well with the eigenstate values. The 
main deviations are near the critical transition values of the energy; from calculations with 

= 100 up to = 10000, we can assert that these deviations should become smaller as 
the number of particles is increased, and should become equal in the thermodynamic limit. 
Hence, we conclude that the stationary states follow the same behavior as the eigenenergy 
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FIG. 6: (Color online) Number of particles in well 1, Ni/N as a function of energy e of the fluid, 
for N = 1000. The frames correspond to A = 0.1, 1.0 and 10.0. Blue solid line corresponds to 
expectation values in energy eigenstates (same as in Fig. [T]). Red line corresponds to the time 
average (mean) of the expectation value in the set of coherent states, Eq.([7]) for ^(e) such that it 
localizes in well 1; purple dotted line corresponds to vr — 6{e) that localizes in well 2. Green dotted 
line corresponds to thermal average values, see Eq.([Io)). 

states and, therefore, that their macroscopic behavior is described by the phase diagram 
shown in Fig. |3j that is, there exists a critical value A^,, such that, for values above such 
a threshold the system exhibits a self-trapping transition as a function of the energy (or 
as a function of the initial state). Additionally, the states may be well characterized by 
both one- and two-body properties; that is, while the one-body properties remain constant 
below the self-trapping transition, two-body variables serve to discriminate among different 
states. Both quantities clearly signal the transition point. Figures [8] and |9] show the standard 
deviation of the values of Ni{t) and C{t)] their smallness justify the characterization of the 
states as stationary. 

The transition from delocalized to self-trapped states exhibits a symmetry-breaking phe- 
nomenon. That is, the Hamiltonian is symmetric under the exchange of the wells, or internal 
states, 1 and 2. However, the stationary states "choose" one of wells to become localized. 
To be precise, one can show that for A > A^, if = and 6 yields a state localized in well 
1, then 9' = TV ~ 9 and = 0, localizes in well 2. This is shown in Fig. |6]for A = 10.0 and 
N = 100, red line signals localization in well 1 and purple line in well 2. 
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FIG. 7: (Color online) Tunneling correlation, C/N"^ as a function of energy e of the fluid, for 
N = 100. The frames correspond to A = 0.1, 1.0 and 10.0. Blue solid line corresponds to 
expectation values in energy eigenstates (same as in Fig. [T]). Red line corresponds to the time 
average (mean) of the expectation value in the set of coherent states, Eq.([7]) for ^(e) such that it 
localizes in well 1; purple dotted line corresponds to vr — 9{e) that localizes in well 2. Green dotted 
line corresponds to thermal average values, see Eq.([Toj). 
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FIG. 8: Standard deviation (fluctuation) of the expectation value the number of particles in well 
1, ANi/N for the set of coherent states, Eq.([7|, as a function of energy. = 1000 and A = 0.1, 
1.0 and 10.0. 

IV. DISCUSSION AND FINAL REMARKS 

We have analyzed one- and two-body properties of the full quantum solution of the 
two-mode Bose-Hubbard fluid. We have discussed the transition from delocalized to self- 
trapped states, which occurs as the energy of the system is increased, provided that the 
pair-interaction strength is above a critical or threshold value. The full quantum solution for 
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FIG. 9: Standard deviation (fluctuation) of the expectation value of the tunneling correlation, 
AC/N"^, for the set of coherent states, Eq.([7|, as a function of energy. N = 1000 and A = 0.1, 1.0 
and 10.0. 

a large number of atoms differs from the mean-field approximation, since the latter predicts 
coherent Josephson-like oscillations while the exact solution shows that those oscillations 
decay to stationary like states in which the system spends most of its time. Since, measuring 
or having access to the A^-body wave function appears as an impractical task in real systems, 
most of our understanding of macroscopic systems is based on knowledge of properties of few 
bodies. In this context, we argue that the decay, or relaxation, to a stationary state can be 
considered as decoherence, even if no interaction with an external environment is included. 
Such a decoherence is a consequence of the interatomic interactions. In a way, this is how 
statistical physics has been developed. That is, one of its basic tenants is the assumption 
that any interacting, macroscopic closed or isolated system, classical or quantum, relaxes to 
an stationary state if left undisturbed [21]. 

Taking the statistical description of the stationary states as an approach to characterize 
them, we have studied the number of particles in well 1 A'^i and the "tunneling correlation" C, 
see Eqs.([5]) and (|6]), as examples of one- and two-body properties. The phase diagram of Fig. 
|3] summarizes our results. For a given values of the number of particles N and the interaction 
strength A, we find a transition from delocalized to self-trapped states as a function of the 
energy e of the fluid, if A > Ac, a critical interaction value. One finds that A^^i remains 
constant, Ni = 0.5N, up to a critical energy value ec(A), where it changes continuously 
up to either a value A^^i = or A^^i = A^, as the energy is further increased. Likewise, the 
transition is registered by C with a cusp at the transition energy, clearly dividing two different 
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types of macroscopic states. This transition appears to us as the reminiscent of a continuous 
phase transition with an spontaneous symmetry breaking mechanism |25j. The behavior of 
A^^i seems analogous to the order parameter of the transition. This is further supported by 
the suggestion, based on our calculations, that while the change of Ni is continuous, its 
derivative appears discontinuous at the transition, and we expect that should become so in 
the thermodynamic limit. The symmetry breaking is also clear: while the Hamiltonian is 
invariant under the interchange "1" and "2" of the wells, stationary states with the same 
energy are not; a set of states localizes in well 1 and another in well 2. This is quite similar 
to the spontaneous magnetization of a ferromagnetic material |24]. 

To conclude, and because of the statistical approach we are putting forward, we would 
like to briefly discuss the role of the temperature in this problem. First, we recall that an 
isolated system is best described by a microcanonical ensemble, and in fact, this is how we 
have described the system. Namely, its macroscopic state is characterized by the number of 
particles and the energy e. Depending on the values of these variables, and if the system is 
in equilibrium, it should have a well defined temperature T = T{N,e). If the ensembles are 
equivalent, one should obtain the same results for the equilibrium states using the canonical 
ensemble with and T given. Due to the simplicity of the system at hand, one can very 
easily calculate the average value of any observable operator A in the canonical ensemble, 

{A) = -J2{<P^\A\<l>n)e-^-/''', (10) 

n=0 

where k is Boltzmann constant and the partition function is 

N 

Z = Y, e"'"/'^- (11) 

n=0 

With this procedure, and assuming that the temperature is a single valued function of 
and e, we can assign a temperature to a stationary state with a given energy e by solving 
the following equation for T, 

1 ^ 

e=-^Y.'ne-'-l'^. (12) 

n=0 



In this way, we can express the thermal averages of A^^i and C by means of Eq.(lO) as a 
function of e and A^, using the temperature T(e, A^) as described above. The green lines in 
Figs. |6] and [7] show one of these examples. The results of the correlation tunneling C seem 
to indicate that there is no agreement between the "microcanonical" and the "canonical" 
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descriptions. In Fig. [g] one finds that the thermal average of Ni is always N/2, whether 
above or below the transition. This is not a discrepancy since the thermal canonical average 
cannot show the spontaneous symmetry breaking; it is analogous to the fact that in the Ising 
model in 2D the average magnetization is always zero: one needs a small external magnetic 
field to break the symmetry and then make it to vanish [TS]. The somewhat perplexing 
problem is in the behavior of the average values of the two-body tunneling correlation C, 
since the peak is not shown in the thermal averages. We do not have an explanation for this 
discrepancy. It may either be that our calculation based on the eigenstates and the coherent 
states do not correspond to a microcanonical ensemble, or that for this system there is simply 
no equivalence between ensembles. The latter option cannot at all be disqualified. In a way 
this system is peculiar since it represents a collection of atoms interacting with the same 
strenght with all the other atoms, provided they are in the same internal state. In other 
words, the fact of considering the two-mode approximation only has given rise to an effective 
Hamiltonian in which the interactions among the particles appear as if they were long range, 
whenever the particles have the same internal state, while the full Hamiltonian represents 



particles with short range interactions. Moreover, it can also be shown, by studying Eq.(12), 
that the energy as a function of temperature and number of particles is not fully extensive, 
namely, it is not of the form e(T, N) = N e(T) with e(T) a function of the energy only, but 
rather, e(T,N) = Ne(T/N). The elucidation of these discrepancies deserves a study on its 
own. 
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